%load alt_calib/eq_mit_mass_shock.mat
clear
%load ../parameter_sandbox_wage_rule/eq_mit_mass.mat 
load ../prod/calibration_v3.5/kimball/eq_mit_entry_cost.mat
%load ../prod/calibration_v3/sandbox/eq_mit_mass

set(0,'defaultaxesfontname','cambria math') % beautify the axes a bit
set(0,'defaultTextFontName', 'cambria math')

Ns = 1 - [0 6.69652130]/100;

adj = sum(eq_mit.mu,1);
scx = pick_scale(adj(2:end), Ns(2));

T = 30;

figure(47)

subplot(3,3,1)
hold on
% plot(scx*(eq_mit.C - 1)); drawnow; title('C'); drawnow;
% xlim([1 T]);
% 
% figure(47)
% 
% subplot(3,3,2)
% hold on
% plot(scx*(eq_mit.D - eq_mit.D(end))); drawnow; title('D'); drawnow;
% xlim([1 T]);
% 
% figure(47)
% 
 Labor       =  eq_mit.l;
 total_labor =  sum(eq_mit.mu.*Labor);
% 
% subplot(3,3,3)
% hold on
% plot((scx*(total_labor - total_labor(end)))/total_labor(end)); title('Total labor demand'); drawnow;
% xlim([1 T]);
% 
% figure(47)
% 
% subplot(3,3,4)
% hold off
% plot(scx*(eq_mit.entry_rate - eq_mit.entry_rate(end)) + eq_mit.entry_rate(end)); title('Entry rate'); drawnow;
% xlim([1 T]);
% 
% figure(47)
% subplot(3,3,5)
% hold off
% plot(sum(eq_mit.exit.*eq_mit.mu)./sum(eq_mit.mu)); title('Exit rate'); drawnow;
% xlim([1 T]);
% 
% figure(47)
% subplot(3,3,6)
% hold off
% plot(mass); title('Total mass'); drawnow;
% xlim([1 T]);
% 
 costs = eq_mit.l;
 prods = repmat(exp(glob.sf(:,2)), 1, options.T);
% a_prods = repmat(sim.A', glob.Nsf, 1);
% prods = prods.*a_prods;
 wages = repmat(eq_mit.W', glob.Nsf, 1);
 mkps  = eq_mit.p./(wages./prods);
 mu = mkps;
% 
 tc    = sum(wages.*costs.*eq_mit.mu);
 num    = sum(wages.*costs.*eq_mit.mu.*mkps);
% 
 mu_cw = num./tc;
% 
% 
% figure(47)
% subplot(3,3,7)
% hold off
% plot(scx*(mu_cw - mu_cw(end)) + mu_cw(end)); title('Cost-weighted markup'); xlim([1 T]); drawnow;
% 
% %
% revs = eq_mit.p.*eq_mit.y;
% 
% tr   = sum(revs.*eq_mit.mu)
% num = sum(revs.*eq_mit.mu.*mkps)/tr;
% 
% mu_sw = num./tr;
% 
% subplot(3,3,8)
% hold off
% plot((mu_sw)); title('Sales-weighted markup'); xlim([1 T]); drawnow;
% 
% 
% print('-dpng', 'figures/irfs_entry.png')


Ys    = repmat(eq_mit.C',glob.Nsf,1);
Z     = sum(eq_mit.y./Ys./prods.*eq_mit.mu);
Z     = Z.^(-1);


%% plot real output, labor demand, and nominal gdp
close all


mass = log_linear_scale(adj, scx);
mass = mass/mass(end);

T = 15;
subplot(2,3,1)
plot(mass, 'LineWidth', 4); xlim([1 T]); ylim([.92 1]); title('Mass of estab.')

subplot(2,3,5)
plot(log_linear_scale(eq_mit.C, scx), 'LineWidth', 4); xlim([1 T]); ylim([.95 1]); title('Output')

subplot(2,3,6)
plot(log_linear_scale(eq_mit.W, scx), 'LineWidth', 4); xlim([1 T]); ylim([.95 1]); title('Wage')

subplot(2,3,4)
plot(log_linear_scale(total_labor, scx)/total_labor(end), 'LineWidth', 4); xlim([1 T]); title('Employment')
ylim([.95 1]);

% subplot(2,3,4)
% ngdp = (eq_mit.entry_rate - eq_mit.entry_rate(end))*scx + eq_mit.entry_rate(end);
% plot(ngdp, 'LineWidth', 4); xlim([1 T]);; ylim([.04 .12]); title('Entry rate')

ngdp = sum(eq_mit.mu.*eq_mit.p.*eq_mit.y);

subplot(2,3,2)
labor_bill = eq_mit.W'.*sum(eq_mit.l.*eq_mit.mu);
labor_share = labor_bill./ngdp;



 plot(log_linear_scale(1./labor_share, scx)/(1./labor_share(end)), 'LineWidth', 4); xlim([1 T]); title('Markup')
%yyaxis right
%mu_cw_s = log_linear_scale(mu_cw, scx);
%plot(mu_cw_s, 'LineWidth', 4);  xlim([1 T]); drawnow;

subplot(2,3,3)
tfp = eq_mit.C'./total_labor;
%tfp = tfp(1:T);
plot(log_linear_scale(tfp, scx)/tfp(end), 'LineWidth', 4); xlim([1 T]);
title('Effective TFP')
ylim([.985 1])

set(gcf,'units','points','position',[10,10,1000,600])
set(findall(gcf,'-property','FontSize'),'FontSize',16)

print('-dpng', 'figures/kimball_entry_cost.png')

